Exploratory Data Analysis

A Demonstration of Tools

Michael Clark https://m-clark.github.io
2020-07-09

Table of Contents


Introduction

In R there are many tools available to help you dive in and explore your data. However, in consulting I still see a lot of people using base R’s table and summary functions, followed by a lot of work to get the result into a more presentable format. My own frustrations led to me creating a package (tidyext) for personal use in this area. While that suits me fine, there are tools that can go much further with little effort. Recently, Staniak & Biecek Staniak and Biecek (2019) wrote an article in the R Journal exploring several of such packages, so I thought I’d try them out for myself, and take others along with me for that ride.

Note that these are first impressions, and I haven’t really dived deeply into any of the packages, and may be missing some key features (apologies to the package authors!). But that is also part of the point for this sort of thing. These aren’t modeling packages, and we have a good idea of what we want in EDA, so these should be easy to pick up and use.

Packages

I’ve updated the article’s table 1, which adds roughly a year’s worth of downloads.

package downloads debut
janitor 495548 2016-10-03
summarytools 199235 2014-08-11
DataExplorer 186456 2016-03-01
visdat 184922 2017-07-11
funModeling 102122 2016-02-07
arsenal 76190 2016-12-30
dlookr 44482 2018-04-27
dataMaid 42799 2017-01-02
inspectdf 23351 2019-04-24
xray 19409 2017-11-22
RtutoR 16272 2016-03-12
ExPanDaR 15464 2018-05-11
exploreR 14459 2016-02-10
SmartEDA 13674 2018-04-06
explore 10125 2019-05-16

A more informative assessment of usage would be in average monthly downloads.

package average_monthly_downloads
janitor 10772.783
visdat 4997.892
DataExplorer 3518.038
summarytools 2767.153
funModeling 1891.148
arsenal 1731.591
dlookr 1588.643
inspectdf 1459.438
dataMaid 995.326
explore 675.000
xray 588.152
ExPanDaR 572.741
SmartEDA 488.357
RtutoR 307.019
exploreR 267.759

Here is a visualization of their growth over time.

Selected packages

Arbitrary criteria

I’ll outline my reasons for selecting some packages to explore and not others. These reasons are somewhat, but not necessarily, arbitrary, and may leave out some viable newer packages. For the data scenario, I am assuming messy data of the sort that might have hundreds of columns of mixed data types, potentially with lots of missingness, attributes that are only applicable to subsets of the data (e.g. branching logic in surveys), etc.1

Here is my criteria for selection:

Things I’m not as concerned about:

In the end, I would like a package that is well put together and will make common tasks easier for me and potentially save me time in creating reports/presentation.

Conceptual organization

Staniak & Biecek note two general phases of data exploration, each with specific tasks, based on the CRISP-DM standard Wirth and Hipp (2000).2

In the first place, I will focus on tools for understanding, particularly description and validity, as they refer to exploration tasks solely as visualization, which is a perk, but something I’m more inclined to do myself.

I’m definitely less interested in the preparation though I will speak about it more towards the end. ‘Cleaning’ in the article refers to mean/median imputation, something I’ve never bothered to do for reasons that have been noted in the statistical literature for a very long time. The other transformations are easy, and probably should be more explicitly documented in your code. Furthermore, in creating derived attributes, things like category merging and standardization depend on the subset of the data used, so it would probably be better to be more explicit than automatic. Also, if things like an automated PCA is viable for your situation, it probably is very simple data (i.e. all variables of the same type), in which case, most of these tools probably won’t add much value to you anyway.

So with that in mind here are the ones I will explore (in alphabetical order):

The Example Dataset

I’ve chosen the heart disease data originally available from the UCI repository (object name hd). It contains a mixture of data types but isn’t too unwieldy, as it’s already been cleaned and has few columns (it’s from my noiris package). For our purposes, I’ve additionally added some random missingness, and created an hd_sample which has only a couple columns to cut down on the display.

Data Description

For data description, we are interested in things like, dimensions of the data, variable types, and maybe even meta-data, like the object’s size in RAM. I also add univariate data summaries as it’s something you’d always want to know and usually report, though these are thought of as data exploration in the article’s delineation.

Preliminary

To give a sense of what my preferences are, consider my own functions. The following actually calls a separate numerical summary function as well as a categorical variable function, and both return ‘tidy’ data frames that can immediately be used for presentation (e.g. via kableExtra) and visualization (e.g. ggplot2), or drilling down to only selections of the output.


library(tidyext)

describe_all(hd, extra = T)  # also describe_all_num, describe_all_cat

$`Numeric Variables`
# A tibble: 8 x 12
  Variable          N   Mean    SD   Min    Q1 Median    Q3   Max `% Missing` Distinct Zeros
  <chr>         <dbl>  <dbl> <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl>       <dbl>    <dbl> <dbl>
1 age             292  54.4   9.15    29   47    55.5  61    77             4       41     0
2 resting_bp      279 132.   17.5     94  120   130   140   200             8       48     0
3 cholesterol     288 246.   51.8    126  211   240   274   564             5      147     0
4 resting_ecg     286   0.51  0.51     0    0     1     1     2             6        3   141
5 max_heartrate   288 150.   23.1     71  133.  154.  167.  202             5       91     0
6 old_peak        284   1.05  1.16     0    0     0.8   1.6   6.2           6       39    91
7 n_vessels       296   0.73  1.02     0    0     0     1     4             2        5   170
8 heart_disease   289   0.53  0.5      0    0     1     1     1             5        2   135

$`Categorical Variables`
# A tibble: 22 x 4
   Variable            Group            Frequency   `%`
   <chr>               <fct>                <int> <dbl>
 1 sex                 male                   199   66 
 2 sex                 female                  95   31 
 3 sex                 <NA>                     9    3 
 4 chest_pain_type     typical angina         143   47 
 5 chest_pain_type     non-anginal pain        82   27 
 6 chest_pain_type     atypical angina         46   15 
 7 chest_pain_type     asymptomatic            19    6 
 8 chest_pain_type     <NA>                    13    4 
 9 fasting_blood_sugar lt_120                 243   80 
10 fasting_blood_sugar gt_120                  42   14.
# … with 12 more rows

I also have options, such as the following.


hd %>% 
  select(age, sex, heart_disease) %>% 
  describe_all(
  digits = 2,
  include_NAcat = FALSE,   # don't include NA as a category
  include_numeric = TRUE,  # allows numeric variables with few levels given in max_levels argument as categorical
  max_levels = 3,
  sort_by_freq = TRUE,
  extra = TRUE
)

$`Numeric Variables`
# A tibble: 2 x 12
  Variable          N  Mean    SD   Min    Q1 Median    Q3   Max `% Missing` Distinct Zeros
  <chr>         <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl>       <dbl>    <dbl> <dbl>
1 age             292 54.4   9.15    29    47   55.5    61    77           4       41     0
2 heart_disease   289  0.53  0.5      0     0    1       1     1           5        2   135

$`Categorical Variables`
# A tibble: 4 x 4
  Variable      Group  Frequency   `%`
  <chr>         <fct>      <int> <dbl>
1 sex           male         199    68
2 sex           female        95    32
3 heart_disease 1            154    53
4 heart_disease 0            135    47

For grouped output, I can use the underlying functions. Here we look at summaries for a couple numerical variables by sex.


hd %>% 
  num_by(
    main_var  = vars(age, cholesterol),
    group_var = sex,
    extra = TRUE
  )

# A tibble: 6 x 13
# Groups:   sex [3]
  sex    Variable        N  Mean    SD   Min    Q1 Median    Q3   Max `% Missing` Distinct Zeros
  <chr>  <chr>       <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl>       <dbl>    <dbl> <dbl>
1 male   age           190  53.5   9      29  46     54     59     77           5       38     0
2 male   cholesterol   190 240.   43.2   126 208    236.   268.   353           5      115     0
3 female age            94  55.4   9.4    34  49.2   56.5   62     76           1       34     0
4 female cholesterol    89 261    65.9   141 214    254    295    564           6       75     0
5 <NA>   age             8  62.1   5.6    54  57     63     67     69          11        6     0
6 <NA>   cholesterol     9 230.   26.9   193 212    231    254    266           0        8     0

The categorical variable functionality is very similar.


hd %>% 
  cat_by(
    main_var  = chest_pain_type,
    group_var = sex
  )

# A tibble: 14 x 5
# Groups:   sex [3]
   sex    chest_pain_type      N `% of Total` `% of sex`
   <chr>  <chr>            <int>        <dbl>      <dbl>
 1 female asymptomatic         3        0.990       3.16
 2 female atypical angina     17        5.61       17.9 
 3 female non-anginal pain    33       10.9        34.7 
 4 female typical angina      39       12.9        41.1 
 5 female <NA>                 3        0.990       3.16
 6 male   asymptomatic        14        4.62        7.04
 7 male   atypical angina     28        9.24       14.1 
 8 male   non-anginal pain    47       15.5        23.6 
 9 male   typical angina     100       33.0        50.3 
10 male   <NA>                10        3.30        5.03
11 <NA>   asymptomatic         2        0.660      22.2 
12 <NA>   atypical angina      1        0.330      11.1 
13 <NA>   non-anginal pain     2        0.660      22.2 
14 <NA>   typical angina       4        1.32       44.4 

As these are tidy tibbles, they are essentially ready for presentation.


describe_all_num(hd) %>% 
  kableExtra::kable()
Variable N Mean SD Min Q1 Median Q3 Max % Missing
age 292 54.36 9.15 29 47.00 55.5 61.00 77.0 4
resting_bp 279 131.91 17.49 94 120.00 130.0 140.00 200.0 8
cholesterol 288 245.95 51.84 126 211.00 240.0 274.00 564.0 5
resting_ecg 286 0.51 0.51 0 0.00 1.0 1.00 2.0 6
max_heartrate 288 149.84 23.11 71 132.75 153.5 167.25 202.0 5
old_peak 284 1.05 1.16 0 0.00 0.8 1.60 6.2 6
n_vessels 296 0.73 1.02 0 0.00 0.0 1.00 4.0 2
heart_disease 289 0.53 0.50 0 0.00 1.0 1.00 1.0 5

These functions serve most of my needs for initial peeking at the data. They return a tibble/data.frame class object that makes for easy presentation and visualization. The underlying code uses the widely used and tested tidyverse packages, and mostly adheres to standard programming conventions. This is the same sort of thing I’m looking for.

arsenal

We begin alphabetically with the arsenal package. Here we use tableby for a generic summary as well as grouped summary. The result is markdown, so for presentation in an R Markdown document, one must use the chunk option results='asis'.


library(arsenal)

hd_sample %>%
  tableby( ~ ., data = .) %>%
  summary()
Overall (N=303)
age
   N-Miss 11
   Mean (SD) 54.356 (9.149)
   Range 29.000 - 77.000
cholesterol
   N-Miss 15
   Mean (SD) 245.955 (51.844)
   Range 126.000 - 564.000
resting_bp
   N-Miss 24
   Mean (SD) 131.907 (17.488)
   Range 94.000 - 200.000
sex
   N-Miss 9
   female 95 (32.3%)
   male 199 (67.7%)
chest_pain_type
   N-Miss 13
   asymptomatic 19 (6.6%)
   atypical angina 46 (15.9%)
   non-anginal pain 82 (28.3%)
   typical angina 143 (49.3%)
heart_disease
   N-Miss 14
   Mean (SD) 0.533 (0.500)
   Range 0.000 - 1.000

hd_sample %>%
  tableby(sex ~ ., data = .) %>%
  summary()
female (N=95) male (N=199) Total (N=294) p value
age 0.092
   N-Miss 1 9 10
   Mean (SD) 55.436 (9.356) 53.495 (8.987) 54.137 (9.140)
   Range 34.000 - 76.000 29.000 - 77.000 29.000 - 77.000
cholesterol 0.001
   N-Miss 6 9 15
   Mean (SD) 261.022 (65.921) 239.658 (43.227) 246.473 (52.397)
   Range 141.000 - 564.000 126.000 - 353.000 126.000 - 564.000
resting_bp 0.325
   N-Miss 8 16 24
   Mean (SD) 133.506 (19.862) 131.240 (16.467) 131.970 (17.627)
   Range 94.000 - 200.000 94.000 - 192.000 94.000 - 200.000
chest_pain_type 0.098
   N-Miss 3 10 13
   asymptomatic 3 (3.3%) 14 (7.4%) 17 (6.0%)
   atypical angina 17 (18.5%) 28 (14.8%) 45 (16.0%)
   non-anginal pain 33 (35.9%) 47 (24.9%) 80 (28.5%)
   typical angina 39 (42.4%) 100 (52.9%) 139 (49.5%)
heart_disease < 0.001
   N-Miss 5 9 14
   Mean (SD) 0.733 (0.445) 0.442 (0.498) 0.536 (0.500)
   Range 0.000 - 1.000 0.000 - 1.000 0.000 - 1.000

Here is an example of categorical-only output using freqlist.


with(hd, table(sex, chest_pain_type)) %>% 
  freqlist() %>% 
  summary()
sex chest_pain_type Freq Cumulative Freq Percent Cumulative Percent
female asymptomatic 3 3 1.07 1.07
atypical angina 17 20 6.05 7.12
non-anginal pain 33 53 11.74 18.86
typical angina 39 92 13.88 32.74
male asymptomatic 14 106 4.98 37.72
atypical angina 28 134 9.96 47.69
non-anginal pain 47 181 16.73 64.41
typical angina 100 281 35.59 100.00

Other options include a function for doing pairwise comparisons for a repeated measure3, comparing datasets, and doing a bunch of bivariate regressions.

Pros

The tableby summary is essentially a ‘Table 1’, which is an unfortunately named display of descriptive stats for a sample of a given study4. My clients often want that, and I’ve actually used a package specifically geared towards providing it just to save me the headache (tableOne), so I definitely find that aspect useful. However, Table 1’s almost invariably have needless statistical output or analysis, and are not a model for reporting that I would choose to go for. As you can see, the layout is not going to be viable for more than a few variables, though it is common practice to present them that way in journal articles regardless of verbosity/legibility.

Issues

I’m not thrilled with markdown output, as I have little control over it, but it’s fine. Also, the default layout is not so succinct, and it appears it’s trying to emulate SAS functions, which are not models for presentation in my opinion. One can use a function to write the output of a single table to html, but I’d rather it just be amenable to the document I’m already creating, or create a report for me of all tables. Lastly, I’m not crazy about using the summary function to get the output. The underlying list objects are not ‘print ready’, so using summary as an argument with default of TRUE would make more sense to me design-wise.

DataExplorer

Now we move to DataExplorer. Let’s introduce the package with the introduce function. I already like this, as it provides good info that extends what you’d get with str or glimpse.


library(DataExplorer)

introduce(hd)

# A tibble: 1 x 9
   rows columns discrete_columns continuous_colu… all_missing_col… total_missing_v… complete_rows total_observati…
  <int>   <int>            <int>            <int>            <int>            <int>         <int>            <int>
1   303      14                6                8                0              214           151             4242
# … with 1 more variable: memory_usage <dbl>

It also can display this information visually in two different ways. I like this, but can’t say I’d ever have a reason to actually use it.


plot_intro(hd)


plot_str(list(hd = hd, gapminder = gapminder_2019))

We can focus on the missingness, which is nice, but I don’t find a use for bar simple bar charts, as they actually make what is only a single row of information harder to parse, and without any visual interest.


plot_missing(hd_sample)

DataExplorer can also plot distributions, e.g. bar and histograms, for the actual values of the categorical/discrete variables.


plot_bar(hd)

But this sort of thing only takes a couple lines of ggplot to do on your own, which you can then customize far more easily. At least DataExplorer can save you the trouble and sorts the result in an ordered fashion by default.


hd %>% 
  select_if(is.character) %>% 
  pivot_longer(everything(), names_to = 'variable', values_to = 'value') %>% 
  ggplot(aes(x = value)) +
  geom_bar() +
  coord_flip() +
  facet_wrap(~ variable, scales = 'free') 

Similarly there are QQ plots, scatterplots and more. I always try to get people to try correlation plots in lieu of large correlation matrices, and DataExplorer provides this. The nice thing is that it will automatically create indicator variables for levels of categorical variables, but beyond that there are issues. For one, it’s diagonal is reversed from the typical presentation of correlation matrices5. If you don’t drop the missing, the plot isn’t as useful, because it will include NA as a factor level, but to its credit DataExplorer has an option to only focus on continuous or discrete output. In addition, ‘centered’ is a rather odd choice for alignment of the x axis in my opinion, so this would take additional work to make presentable.


plot_correlation(
  hd,
  ggtheme = theme_minimal(),
  cor_args = list("use" = "pairwise.complete.obs")
)

There are packages that do this specifically, such as heatmaply. Before these came around I already had my own function, that also does an internal factor analysis to sort the variables, and produces an interactive result. So while so this aspect of DataExplorer might be useful to others, it doesn’t appeal much to me, especially given the other issues.


hd %>% 
  select_if(is.numeric) %>% 
  cor(use = 'pair') %>% 
  visibly::corr_heat()

If the various options of output, or only some of them, appeal to you, they can all be nicely wrapped up in an automatic report. Various settings for each type of output can be set with an additional function (configure_report) or just passing a list of arguments, including which functions to use, ggplot2 theme, etc. You can see the report here, but I show a screenshot.


create_report(
  hd,
  y = 'heart_disease',
  output_dir   = 'other_docs',
  output_file  = 'data_explorer_report.html',
  report_title = 'My Data Description'
)

Pros

I think many would like at least some functionality in DataExplorer, as well as many of the visualizations. The ease with which to generate a report should also be sufficient for anyone’s personal use, and with some tweaking, presentation to others. It also uses data.table under the hood, so likely can handle large data with efficiency.

Not covered here, but DataExplorer also has functionality for feature processing and engineering, for example, collapsing sparse categories, dummy coding, etc.

Issues

The issues I have are pretty minor with this one aside from unnecessary statistical analysis and visualization choices. I would also recommend pryr::object_size rather than base R’s function.

SmartEDA

I wanted to look at the SmartEDA package because the figure in the article was of a clean report. Let’s start our exploration with the basic ExpData function.


library(SmartEDA)

ExpData(hd)

                                Descriptions       Obs
1                         Sample size (Nrow)       303
2                    No. of Variables (Ncol)        14
3                   No. of Numeric Variables         8
4                    No. of Factor Variables         0
5                      No. of Text Variables         6
6                   No. of Logical Variables         0
7                    No. of Unique Variables         0
8                      No. of Date Variables         0
9   No. of Zero variance Variables (Uniform)         0
10     %. of Variables having complete cases    0% (0)
11 %. of Variables having <50% missing cases 100% (14)
12 %. of Variables having >50% missing cases    0% (0)
13 %. of Variables having >90% missing cases    0% (0)

We can look at ExpNumStat to get some basic stats for numeric variables.


ExpNumStat(hd, round = 1)

          Vname Group  TN nNeg nZero nPos NegInf PosInf NA_Value Per_of_Missing     sum min   max  mean median   SD  CV
1           age   All 303    0     0  292      0      0       11            3.6 15872.0  29  77.0  54.4   55.5  9.1 0.2
3   cholesterol   All 303    0     0  288      0      0       15            5.0 70835.0 126 564.0 246.0  240.0 51.8 0.2
4 max_heartrate   All 303    0     0  288      0      0       15            5.0 43155.0  71 202.0 149.8  153.5 23.1 0.2
5      old_peak   All 303    0    91  193      0      0       19            6.3   296.8   0   6.2   1.0    0.8  1.2 1.1
2    resting_bp   All 303    0     0  279      0      0       24            7.9 36802.0  94 200.0 131.9  130.0 17.5 0.1
   IQR Skewness Kurtosis
1 14.0     -0.2     -0.6
3 63.0      1.2      4.7
4 34.5     -0.5     -0.1
5  1.6      1.3      1.6
2 20.0      0.7      0.9

We can look at ExpNumStat to get some basic stats for grouped output also, and set various options.


ExpNumStat(
  hd_sample,
  by = "GA",
  gp = "sex",
  Qnt = c(.1, .9),
  Outlier = TRUE,
  round = 1
)

         Vname      Group  TN nNeg nZero nPos NegInf PosInf NA_Value Per_of_Missing   sum min  max  mean median   SD
1          age    sex:All 303    0     0  292      0      0       11            3.6 15872  29   77  54.4   55.5  9.1
4          age   sex:male 199    0     0  190      0      0        9            4.5 10164  29   77  53.5   54.0  9.0
7          age sex:female  95    0     0   94      0      0        1            1.1  5211  34   76  55.4   56.5  9.4
10         age     sex:NA   0    0     0    0      0      0        0            NaN     0 Inf -Inf   NaN     NA   NA
2  cholesterol    sex:All 303    0     0  288      0      0       15            5.0 70835 126  564 246.0  240.0 51.8
5  cholesterol   sex:male 199    0     0  190      0      0        9            4.5 45535 126  353 239.7  235.5 43.2
8  cholesterol sex:female  95    0     0   89      0      0        6            6.3 23231 141  564 261.0  254.0 65.9
11 cholesterol     sex:NA   0    0     0    0      0      0        0            NaN     0 Inf -Inf   NaN     NA   NA
3   resting_bp    sex:All 303    0     0  279      0      0       24            7.9 36802  94  200 131.9  130.0 17.5
6   resting_bp   sex:male 199    0     0  183      0      0       16            8.0 24017  94  192 131.2  130.0 16.5
9   resting_bp sex:female  95    0     0   87      0      0        8            8.4 11615  94  200 133.5  132.0 19.9
12  resting_bp     sex:NA   0    0     0    0      0      0        0            NaN     0 Inf -Inf   NaN     NA   NA
    CV  IQR Skewness Kurtosis   10%   90% LB.25% UB.75% nOutliers
1  0.2 14.0     -0.2     -0.6  42.0  66.0   26.0   82.0         0
4  0.2 13.0     -0.2     -0.5  41.9  65.0   26.5   78.5         0
7  0.2 12.8     -0.3     -0.6  42.0  66.0   30.1   81.1         0
10  NA   NA      NaN      NaN    NA    NA     NA     NA         0
2  0.2 63.0      1.2      4.7 188.0 308.3  116.5  368.5         5
5  0.2 60.5      0.1     -0.3 184.9 299.1  117.2  359.2         0
8  0.3 81.0      1.4      4.0 197.0 327.6   92.5  416.5         2
11  NA   NA      NaN      NaN    NA    NA     NA     NA         0
3  0.1 20.0      0.7      0.9 110.0 152.4   90.0  170.0         8
6  0.1 20.0      0.7      0.7 110.0 150.0   90.0  170.0         3
9  0.1 20.0      0.7      0.8 109.2 160.0   90.0  170.0         5
12  NA   NA      NaN      NaN    NA    NA     NA     NA         0

There is also some visualization of relationships with a given target variable.


ExpNumViz(hd_sample, target = 'cholesterol')

[[1]]


[[2]]

Here are the default categorical data summaries. This is a nice and clean data.frame presentation.


ExpCTable(hd_sample) 

          Variable            Valid Frequency Percent CumPercent
1              sex           female        95   31.35      31.35
2              sex             male       199   65.68      97.03
3              sex               NA         9    2.97     100.00
4              sex            TOTAL       303      NA         NA
5  chest_pain_type     asymptomatic        19    6.27       6.27
6  chest_pain_type  atypical angina        46   15.18      21.45
7  chest_pain_type               NA        13    4.29      25.74
8  chest_pain_type non-anginal pain        82   27.06      52.80
9  chest_pain_type   typical angina       143   47.19      99.99
10 chest_pain_type            TOTAL       303      NA         NA
11   heart_disease                0       135   44.55      44.55
12   heart_disease                1       154   50.83      95.38
13   heart_disease               NA        14    4.62     100.00
14   heart_disease            TOTAL       303      NA         NA

As there was with numeric variables, there is also visualization for the categorical variables.


ExpCatViz(hd_sample) 

[[1]]


[[2]]


[[3]]

As far as reporting, SmartEDA also provides this functionality. There are some options you can fiddle with, like using your own template, changing the default theme, etc.


ExpReport(
  hd,
  theme   = visibly::theme_clean(),
  op_dir  = 'other_docs/',
  op_file = 'smarteda.html'
)

Pros

SmartEDA is fairly intuitive to use. It returns a data frame, and some of the less verbose output is quite ready to go. It can also generate a report in automatic fashion.

Issues

A lot of this isn’t very useful to me, such as parallel coordinate plots, ‘outlier’ analysis, etc. I’m not crazy about the naming conventions, both for the functions and arguments (every single function in the package begins with Exp). The default color schemes for some output is simply not viable for presentation, and the report doesn’t really have any options for controlling the output like DataExplorer did. I altos had issues with attempting to get it to work outside of the current working directory for the initial document creation. And changing the sn or sc options, which I’m not sure why I would only want a sample of plots rather than specifying which plots I’d want specifically, didn’t appear to actually change the resulting document.

summarytools

The summarytools package provides four main functions to work with, but I’m going to skip those and go straight to the tool that uses those key functions and puts their results into a very nice presentation.


library(summarytools)

dfSummary(
  hd,
  varnumbers = FALSE,
  round.digits = 2,
  plain.ascii = FALSE,
  style = "grid",
  graph.magnif = .33,
  valid.col = FALSE,
  tmp.img.dir = "img"
)

Data Frame Summary

hd

Dimensions: 303 x 14
Duplicates: 0

Variable Stats / Values Freqs (% of Valid) Graph Missing
age
[numeric]
Mean (sd) : 54.4 (9.1)
min < med < max:
29 < 55.5 < 77
IQR (CV) : 14 (0.2)
41 distinct values 11
(3.63%)
sex
[character]
1. female
2. male
95 (32.3%)
199 (67.7%)
9
(2.97%)
chest_pain_type
[character]
1. asymptomatic
2. atypical angina
3. non-anginal pain
4. typical angina
19 ( 6.6%)
46 (15.9%)
82 (28.3%)
143 (49.3%)
13
(4.29%)
resting_bp
[numeric]
Mean (sd) : 131.9 (17.5)
min < med < max:
94 < 130 < 200
IQR (CV) : 20 (0.1)
48 distinct values 24
(7.92%)
cholesterol
[numeric]
Mean (sd) : 246 (51.8)
min < med < max:
126 < 240 < 564
IQR (CV) : 63 (0.2)
147 distinct values 15
(4.95%)
fasting_blood_sugar
[character]
1. gt_120
2. lt_120
42 (14.7%)
243 (85.3%)
18
(5.94%)
resting_ecg
[numeric]
Mean (sd) : 0.5 (0.5)
min < med < max:
0 < 1 < 2
IQR (CV) : 1 (1)
0 : 141 (49.3%)
1 : 143 (50.0%)
2 : 2 ( 0.7%)
17
(5.61%)
max_heartrate
[numeric]
Mean (sd) : 149.8 (23.1)
min < med < max:
71 < 153.5 < 202
IQR (CV) : 34.5 (0.2)
91 distinct values 15
(4.95%)
exer_angina
[character]
1. no
2. yes
193 (67.0%)
95 (33.0%)
15
(4.95%)
old_peak
[numeric]
Mean (sd) : 1 (1.2)
min < med < max:
0 < 0.8 < 6.2
IQR (CV) : 1.6 (1.1)
39 distinct values 19
(6.27%)
slope
[character]
1. flat
2. negative
3. positive
130 (45.6%)
134 (47.0%)
21 ( 7.4%)
18
(5.94%)
n_vessels
[numeric]
Mean (sd) : 0.7 (1)
min < med < max:
0 < 0 < 4
IQR (CV) : 1 (1.4)
0 : 170 (57.4%)
1 : 65 (22.0%)
2 : 36 (12.2%)
3 : 20 ( 6.8%)
4 : 5 ( 1.7%)
7
(2.31%)
defect
[character]
1. fixed_defect
2. normal
3. reversible_defect
153 (53.9%)
18 ( 6.3%)
113 (39.8%)
19
(6.27%)
heart_disease
[numeric]
Min : 0
Mean : 0.5
Max : 1
0 : 135 (46.7%)
1 : 154 (53.3%)
14
(4.62%)

This is exactly what I want- basic, not overwhelming and redundant information, a usable data frame object, simple visuals to enhance the output without adding to the total information you have to parse (and which can be turned off), and basic categorical information. In one function, I have pretty much all I’d need, but with control to tweak as necessary.


hd_sample %>% 
  group_by(sex) %>% 
  dfSummary(
    varnumbers = FALSE,
    round.digits = 2,
    plain.ascii = FALSE,
    na.col = FALSE,
    style = "grid",
    graph.magnif = .33,
    valid.col = FALSE,
    tmp.img.dir = "/tmp"
  )

Data Frame Summary

hd_sample

Group: sex = female
Dimensions: 95 x 6
Duplicates: 0

Variable Stats / Values Freqs (% of Valid) Graph
age
[numeric]
Mean (sd) : 55.4 (9.4)
min < med < max:
34 < 56.5 < 76
IQR (CV) : 12.8 (0.2)
34 distinct values
cholesterol
[numeric]
Mean (sd) : 261 (65.9)
min < med < max:
141 < 254 < 564
IQR (CV) : 81 (0.3)
75 distinct values
resting_bp
[numeric]
Mean (sd) : 133.5 (19.9)
min < med < max:
94 < 132 < 200
IQR (CV) : 20 (0.1)
34 distinct values
sex
[character]
1. female 95 (100.0%)
chest_pain_type
[character]
1. asymptomatic
2. atypical angina
3. non-anginal pain
4. typical angina
3 ( 3.3%)
17 (18.5%)
33 (35.9%)
39 (42.4%)
heart_disease
[numeric]
Min : 0
Mean : 0.7
Max : 1
0 : 24 (26.7%)
1 : 66 (73.3%)

Group: sex = male
Dimensions: 199 x 6
Duplicates: 0

Variable Stats / Values Freqs (% of Valid) Graph
age
[numeric]
Mean (sd) : 53.5 (9)
min < med < max:
29 < 54 < 77
IQR (CV) : 13 (0.2)
38 distinct values
cholesterol
[numeric]
Mean (sd) : 239.7 (43.2)
min < med < max:
126 < 235.5 < 353
IQR (CV) : 60.5 (0.2)
115 distinct values
resting_bp
[numeric]
Mean (sd) : 131.2 (16.5)
min < med < max:
94 < 130 < 192
IQR (CV) : 20 (0.1)
40 distinct values
sex
[character]
1. male 199 (100.0%)
chest_pain_type
[character]
1. asymptomatic
2. atypical angina
3. non-anginal pain
4. typical angina
14 ( 7.4%)
28 (14.8%)
47 (24.9%)
100 (52.9%)
heart_disease
[numeric]
Min : 0
Mean : 0.4
Max : 1
0 : 106 (55.8%)
1 : 84 (44.2%)

Group: sex = NA
Dimensions: 9 x 6
Duplicates: 0

Variable Stats / Values Freqs (% of Valid) Graph
age
[numeric]
Mean (sd) : 62.1 (5.6)
min < med < max:
54 < 63 < 69
IQR (CV) : 10 (0.1)
54 : 1 (12.5%)
57 : 2 (25.0%)
62 : 1 (12.5%)
64 : 1 (12.5%)
67 : 2 (25.0%)
69 : 1 (12.5%)
cholesterol
[numeric]
Mean (sd) : 229.9 (26.9)
min < med < max:
193 < 231 < 266
IQR (CV) : 42 (0.1)
193 : 1 (11.1%)
201 : 1 (11.1%)
212 : 2 (22.2%)
231 : 1 (11.1%)
239 : 1 (11.1%)
254 : 1 (11.1%)
261 : 1 (11.1%)
266 : 1 (11.1%)
resting_bp
[numeric]
Mean (sd) : 130 (13.3)
min < med < max:
110 < 125 < 152
IQR (CV) : 16 (0.1)
110 : 1 (11.1%)
120 : 1 (11.1%)
124 : 2 (22.2%)
125 : 1 (11.1%)
130 : 1 (11.1%)
140 : 1 (11.1%)
145 : 1 (11.1%)
152 : 1 (11.1%)
sex
[character]
All NA’s
chest_pain_type
[character]
1. asymptomatic
2. atypical angina
3. non-anginal pain
4. typical angina
2 (22.2%)
1 (11.1%)
2 (22.2%)
4 (44.4%)
heart_disease
[numeric]
Min : 0
Mean : 0.4
Max : 1
0 : 5 (55.6%)
1 : 4 (44.4%)

Pros

The summarytools package provides a single function that produces a ready-to-present table within whatever document I’m already using. Very nice!

Issues

The only nitpicky stuff I have with this package is that the ctable function isn’t useful to me, and some of the numeric description is not really necessary. I had issues with the trying to use the suggested tmp directory for the image files, but found it easier to just use a project folder.

Data Validity

As a reminder, data validity is more about providing checks on the data rather than summarizing it per se. These are packages more geared toward spotting or dealing with data issues.

dataMaid

The primary utility of dataMaid is examination of data consistency, but it will also provide basic summaries too. In the end it creates a useful report as well. I came across this via the author’s presentation at the Ann Arbor R User Group, and have actually used it before for a project with success. The result can be found here, though here is a screenshot.


library(dataMaid)

dataMaid::makeDataReport(
  hd,
  output = 'html',
  file = 'other_docs/dataMaid_report',
  replace = TRUE,
  maxDecimals = 1
) 

As you can see, it provides an overall summary of the data frame including variable types, unique values, missing percentage, and potential problems. Problems mostly regard outliers, which is fine to inspect, but arbitrarily set.

Pros

dataMaid is highly customizable, but the default is already a great way to get a good sense of whether your data is doing what it’s supposed to do.

Issues

I initially couldn’t get anything but a pdf document even though I clearly specified html as per the documentation, so I’m not sure what’s going on there, but I eventually sorted it out. The default outlier check, while better than most would do, seemed too sensitive, but this is nitpicky. In the past I’ve found it be slow in rendering the output for a larger data set, but this isn’t something you would usually have to run more than once after getting your options squared away.

janitor

Sadly, I still regularly receive data from Excel (and SPSS), which means I have to worry about things that I shouldn’t have to worry about, like whether dates will actually be treated appropriately, if the column names are usable, or whether the imported file contains 10000 empty rows because someone accidentally hit the space bar.

The janitor package provides a few simple tools that many will probably find useful at some point, and I regularly use its remove_empty function after importing anything from Excel. As an example I’ll add an empty column with an Excel-like name, a duplicated row, and a useless column with a constant value (all very typical in Excel files I receive).


library(janitor)

hd2 = hd

hd2$`test this - ridiculous (name)` = NA           # a terribly formatted name that has no values
hd2$`why.is-this_here` = 'same value everywhere'   # constant column
hd2[nrow(hd2) + 1,] = hd2[1,]                      # add a duplicated row

hd2 %>% 
  clean_names() %>% 
  colnames()

 [1] "age"                       "sex"                       "chest_pain_type"           "resting_bp"               
 [5] "cholesterol"               "fasting_blood_sugar"       "resting_ecg"               "max_heartrate"            
 [9] "exer_angina"               "old_peak"                  "slope"                     "n_vessels"                
[13] "defect"                    "heart_disease"             "test_this_ridiculous_name" "why_is_this_here"         

And here is the remove_* functionality.


hd3 = hd2 %>% 
  remove_empty() %>% 
  remove_constant() 

# useless columns/rows removed 
colnames(hd3)

 [1] "age"                 "sex"                 "chest_pain_type"     "resting_bp"          "cholesterol"        
 [6] "fasting_blood_sugar" "resting_ecg"         "max_heartrate"       "exer_angina"         "old_peak"           
[11] "slope"               "n_vessels"           "defect"              "heart_disease"      

c(nrow(hd2), nrow(hd3))

[1] 304 304

Just a couple little things like that are very useful. Otherwise, there is functionality for dates, tables, etc. that some might find useful in a pinch, but one would probably will have more with the other summary functions in packages we’ve demoed, and more advanced users may simply just use lubridate, stringr and similar packages directly.

Pros

Provides useful functionality not seen in the other packages.

Issues

Not a whole lot going on here relative to some of the other packages, but what is there is likely to be useful to many.

visdat

The package visdat6 is, as its name implies, purely for visualization, and this includes missingness, correlation, and more. We can start by visualizing the missing data.


library(visdat)

vis_dat(hd_sample)

We can look at variable types along with the missingness. I like the integer/double mix, which might point to an issue if all the data should be integer valued.


vis_guess(hd)

We can visualize the correlations as we did with DataExplorer, and default is pairwise. This is one of the cleaner correlation plots I’ve come across, but there is no way to order it meaningfully.


hd %>% 
  select_if(is.numeric) %>% 
  vis_cor()

One nice feature I haven’t seen elsewhere is to visualize a given expected value across a data set.


hd %>% 
  select(sex) %>% 
  vis_expect(expectation =  ~ .x == 'male')

We can also compare whole data sets. Here the differences are with missing values in one and not the other7.


hd2 = noiris::heart_disease %>% 
  select(colnames(hd_sample)) %>% 
  mutate_if(is.factor, as.character)

vis_compare(
  hd_sample,
  hd2
)

Pros

This package adds some useful functionality we haven’t seen. In addition, the underlying code adheres to open science standards. There is a complementary package naniar for more ways to deal with missing data.

Issues

This isn’t for numeric description, so can only supplement typical EDA as we’ve had before. You may have to do some pre-processing for some functions unlike with other tools (e.g. subsetting to numeric). But these are minor issues at best given the package’s purpose.

Data Cleaning and Transformations

As mentioned previously, I’m not interested in using these packages for doing this, and based on Staniak and Biecek, only a couple of them provide this, DataExplorer being the only package that we’ve explored here. Mean/median/mode imputation is a great way to attenuate correlation, so I’m not interested in doing that. I don’t even know what an ‘outlier’ is outside of a modeling framework, so I’m definitely not interested in doing something about extreme values based on a univariate analysis8. Discretizing numeric variables should almost never be done9, so that functionality is not desirable. And between base R functions like scale and log, or the rescale function from scales, lumping factors and similar via forcats, I don’t really need another package for this sort of thing.

Summary

Staniak & Biecek Staniak and Biecek (2019) provide a nice summary of many packages that can be used for automated exploratory data analysis. I wanted to explore and demonstrate some of these in more detail. My impression is that many of these packages have notably useful functionality, though a lot of it may be overlapping, or superfluous from a more serious analytical standpoint. In addition, I personally wouldn’t use the word automated to describe most of the functionality for most of these packages, as multiple specific functions, possibly with different options for each, would need to be used to ultimately generate a data driven product. That’s not necessarily a bad thing, as one should know what’s in the data if you want to make a decision based on it.

Here is a rough list of features I made based on quick inspections.

Table 1: See Table 4 in Staniak & Biecek for more details.
Package Ready to use output Code Quality Visualization Report Generation Dedicated Website
arsenal 👍
DataExplorer 👍 👍 👍 👍 👍
SmartEDA 👍 👍 👍 👍
summarytools 👍 👍
visdat 👍 👍 👍 👍
dataMaid 👍 👍
janitor 👍 👍 👍

In general I would imagine I’d use summarytools or DataExplorer, with visdat and janitor to fill in additional information or help with the data processing. Interestingly, these also appear to be the most popular packages under consideration10, so it’s good to know they are more widely used. Any of them might be useful for your specific needs, particularly if you’re willing to invest in some of the more customizable ones, so take a spin with any you like.

Description packages

arsenal

DataExplorer

SmartEDA

summarytools

Validation packages

visdat

dataMaid

janitor

Exercises

Find a data set you like and do the following. I suggest using a manageable one or some subset so that you don’t get an overwhelming amount of output.

Staniak, Mateusz, and Przemyslaw Biecek. 2019. “The Landscape of R Packages for Automated Exploratory Data Analysis.” The R Journal. https://journal.r-project.org/archive/2019/RJ-2019-033/index.html.

Wirth, Rüdiger, and Jochen Hipp. 2000. “CRISP-Dm: Towards a Standard Process Model for Data Mining.” In Proceedings of the 4th International Conference on the Practical Applications of Knowledge Discovery and Data Mining, 29–39. Springer-Verlag London, UK.


  1. However, for our demonstration I’ll be using something a little more wieldy.↩︎

  2. To be honest, I wasn’t familiar with the CRoss Industry Standard Process for Data Mining until reading the article citing it. I don’t get the impression any particular methodology is actually consciously thought about by the vast majority of practicing data scientists, but it’s useful in providing a framework for the content here.↩︎

  3. The data we’re using isn’t paired and their documentation doesn’t provide a working example.↩︎

  4. I see ‘Table 1’ used mostly by folks in medical fields, who subsequently place it as table 2, 3 or whatever as would normally be the case.↩︎

  5. I suspect this may have more to do with ggplot than what is desired. But I also found that the orientation/alignment of names was different for the report versus the inline presentation.↩︎

  6. The article by Staniak & Biecek places visdat with description rather than validity, but the data frame comparison, expectation investigation, and missing exploration seem more validity issues than data summaries.↩︎

  7. The arsenal package also provided this functionality. It’s quite verbose, so I didn’t show it, but if you need a more in depth comparison it can be recommended.↩︎

  8. Outliers are an indicator of model weakness/limitation/failure, not a data problem.↩︎

  9. The only times I do this is with an already discrete variable, for example, going from 7 values to 5 values, or coarsening that might be applied in very large data situations for computational efficiency. Issues with discretizing are noted here, and while there are better ways to do it if you feel you must, you probably won’t have seen it done in practice, and the chosen method likely wouldn’t hold across repeated data.↩︎

  10. The janitor package, which was not explored in the original article, is actually notably more popular than any of the other packages.↩︎